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ABSTRACT 

Cosniological simulations of structure formation follow the coUisionless evolution of 
dark matter starting from a nearly homogeneous field at early times down to the highly 
clustered configuration at redshift zero. The density field is sampled by a number of 
particles in number infinitely smaller than those believed to be its actual components 
and this limits the mass and spatial scales over which we can trust the results of a 
simulation. Softening of the gravitational force is introduced in coUisionless simula- 
tions to limit the importance of close encounters between these particles. The scale 
of softening is generally fixed and chosen as a compromise between the need for high 
spatial resolution and the need to limit the particle noise. In the scenario of cosmo- 
logical simulations, where the density field evolves to a highly inhomogeneous state, 
this compromise results in an appropriate choice only for a certain class of objects, 
the others being subject to either a biased or a noisy dynamical description. We have 
implemented adaptive gravitational softening lengths in the cosmological simulation 
code GADGET-3; the formalism allows the softening scale to vary in space and time 
according to the density of the environment, at the price of modifying the equation 
of motion for the particles in order to be consistent with the new dependencies intro- 
duced in the system's Lagrangian. We have applied the technique to a number of test 
cases and to a set of cosmological simulations of structure formation. We conclude 
that the use of adaptive softening enhances the clustering of particles at small scales, 
a result visible in the amplitude of the correlation function and in the inner profile of 
massive objects, thereby anticipating the results expected from much higher resolution 
simulations. 
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1 INTRODUCTION 

In coUisionless N-body simulations the matter field is rep- 
resented by a number of point-particles in number consid- 
erably reduced with respect to the actual building blocks 
of the system. These particles have no physical counterpart 
and should be simply regarded as Monte-Carlo samplings 
of the probability-density distribution in position and ve- 
locity. When during the evolution two of these particles are 
found at small spatial separations, the gravitational attrac- 
tion between them can become arbitrarily large and follow- 
ing the encounter properly can turn out to be prohibitively 
expensive from the computational point of view; plus, this 
collisional behaviour is an artifact of the granularity of the 
mass distribution in the simulation and would not manifest 
in a intrinsically-smooth, coUisionless system. A way around 
this problem, namely the discreetness of the particle repre- 



E-mail: iannuzzi@mpa-garching.mpg.de 



methods: numerical - cosmology: theory, large-scale struc- 



sentation, is obtained by smoothing the density distribution 
from a collection of spikes to a continuous field; to this pur- 
pose, every particle is assigned a finite volume over which its 
mass is spread. This translates into the gravitational interac- 
tion between the computational particles being "softened" 
at small separations, i.e. being attenuated and prevented 
to diverge. In this way close encounters are impeded and 
the simulation can proceed at a regular pace. The gravita- 
tional smoothing does not considerably affect the artificial 
two-body relaxation though, as this is driven by both close 
and distant encounter s and the former cannot be avoided by 
softening techniques (|Chandrasekhail 19421; ISpitzer fc Hard 
19711; iHernauist fc Barnes! 1 1990l ; lTheislll998l ; iDehnenllJOOll ; 
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Diemand et al.ll2004 ). Suppressing relaxation is primarily a 



condition on the number of sampling pa rticles and depend s 
only weakly on the value of the softening (jPower et al.ll2003l ) . 
The details of this smoothing procedure can be reduced to 
the choice of a softening kernel W{r,h), which determines 
the functional form of the modified density profile of the par- 
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tide, and of a softening length h, which controls the spatial 
extent within which the modification applies. In the classic 
case of "Plummer" softening, the gravitational potential of 
each particle is that of a Plummer sphere whose scale length 
is given by the value of the softening h; this results in the 
following form for the gravitational force F{r) between two 
particles of masses rrii and rrij separated by a distance r: 



F{r) = -G 



(r2 + /l2)3/2 



(1) 



An evident shortcoming of this choice is that the gravita- 
tional force is modified at all sep arations; several works (e.g. 
iDver fc id Il993l : iDehnenI l200lh have shown that a better 
choice is to use kernels with a compact support, meaning 
that the interaction recovers its Newtonian, original form at 
separations grater than the softenin g length. The most com- 
monly used one is the cubic spline of lMonaghan fc Lattanzid 
l|l985l ): 



W{r,h) = -——{ 2(l-q)^ 0.5 <g< 1 (2) 
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where q = r/h, r being the distance from the centre of the 
kernel, i.e. from the particle's position. The particles assume 
the density profile given by Eq. [2] and the gravitational po- 
tential and force fields are modified accordingly. 
It is evident that the introduction of softening leads to an 
unphysical modification of the gravitational interaction at 
interparticle separation below h\ the inverse square law is 
replaced by a gentler interaction, whose strength approaches 
zero in the limit r ^ 0. This results in a systematic misrep- 
resentation of the force at small scales, an effe ct commonly 
referred to as bias lJMerrittlll996l : lDehnenll200ll ). Ideally one 
would like to limit this effect to the smallest possible scales 
by reducing h accordingly; at the same time the softening 
lengths cannot be made arbitrarily small without running 
into the discreetness problem previously addressed. It is then 
clear how the choice of h assumes a critical importance for 
the reliabili ty of the force computation and indeed a number 
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have been devoted to the subject. 
In general there is an overall agreement that no such a thing 
as an "optimal" softening length exists when one has to 
deal with highly inhomogeneous systems, as too big a value 
would degrade the description of over-dense regions whilst 
too small a value would enhance collisionality in under-dense 
environments. Generally, the gravitational softening scale is 
set to a fixed value, indeed decided as a compromise between 
the desired spatial resolution and the need to moderate the 
noise arising from particle-particle interactions. As a matter 
of fact, in most cases this means fine-tuning the softening 
parameter to follow properly the dynamics of a certain class 
of objects - the densest - at a certain time; in a cosmological 
simulations, where the initially uniform matter density field 
evolves to form highly differentiated structures, this implies 
following poorly not only the moderately-dense regions, but 
also the progenitors of the "target", densest objects. 
Ideally it would be desirable to vary the softening scale in 
space and time according to the density evolution of the dif- 
ferent environments: this would remove the aforementioned 
trade-off between spatial resolution and particle noise, thus 
allowing to increase the former and reduce the latter depend- 



ing on the local features of the density field. The problem 
with a varying softening length though, as will be shown in 
Section [21 is that it introduces an additional dependence on 
the spatial position r in the definition of the gravitational 
potential; if this is not accounted for in the derivation of the 
equation of motion, the energetics of the system can manifest 
unwan ted behaviours. In this context, IPrice fc Monaghad 
((2003) (hereafter PM07) have proposed a formalism to adapt 
the gravitational softening lengths in a simulation while re- 
taining conservation of both momentum and energy, which, 
as just mentioned, would be lost when adopting an adaptive 
scheme. In their method the softening lengths are computed 
in the same fashion as for the s moothing length in smoothed 
particle hy drodynamics (SPH - iGingold fc Monaghadll977l : 
ILucvI 19771 1 , where h varies with the local density p according 
to ft oc p~^'^ ■ The adoption of this scheme involves a num- 
ber of modifications in the overall structure of a typical cos- 
mological simulation code, most noticeably the introduction 
of an additional term in the equations of motion. The tech- 
nique has alread y been employed in simul ation codes like the 
TreePM code of iBagla fc Khandail (|2009|) (he r eafter BK09) 
and the TreeSPH EvoL code of iMerlin et al] (|201G| '). BK09 
apply the method to coUisionless simulations of structure 
formation and find an increase in clustering in over-dense 
regions accompanied by a somewhat loss in resolution in 
the under-dense ones; iMerlin et al.l (|201Gr ) performed stan- 
dard hydrodynamical tests in order to test the overall per- 
formance of their code and found good energy- conservation 
properties along with an improved behaviour i n tests like 
the isothermal collapse of lBate fc BurkertI (|l997l ). Although 
we are here concerned about particle-based codes, where the 
gravitational interaction is computed at the particle's loca- 
tion, we not e that grid-based codes employing adaptive mesh 



refinement (IBrvan fc Norman 
Abel et all I2OO0I; iKnebe et al 



1997; Truelove et al 



2001 : Kravtsov et al 



1998 



2002 



Tevssiei) |2002| ; lO'Shea et al.ll2004l ; iQuilisI |2004| ) are intrin 
sically spatially adaptive; these codes perform a recursive 
refinement of an initial regular grid, thereby increasing their 
resolution in regions of interest. For similar reasons as those 
outlined above, conservation of energy is not strictly guar- 
anteed in these codes, as it would not be when using adap- 
tive gravitational softening without changing the equation 
of motion accordingly. 

In this paper we discuss the implementation of the adap- 
tive softening length form alism by PM07 in the cosmological 
simulation co de GADGET jSpringel. Yoshida fc Whit JI2OOII : 
ISpringel2005l ) and the main effects that it has on dark mat- 
ter simulations of the large-scale structure of the Universe. 
Our results qualitatively confirm those of BK09 in over- 
dense regions, but show no drawbacks in under-dense ones; 
in particular, we do not register a significant loss of low-mass 
structures and we considerably improve the results in term 
of particle clustering when comparing to standard simula- 
tions employing fixed softening. 

In Sections [2] and [3] we briefiy describe the technique and 
its implementation in GADGET ; in SectionUwe show two of 
the tests performed to check the implementation and its re- 
sult in controlled scenarios for which we know the solution, 
whereas in Sections O and [S] we apply the technique to cos- 
mological simulations of structure formation. We summarise 
the results and conclude in Section [T] 
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2 THE FORMALISM 

In this section we will briefly report the gist of the algorithm. 
We refer the reader to the pivotal work PM07 for a thorough 
treatment of the problem and for a detailed derivation of the 
new equation of motion. 

The Lagrangian for a system of A'' particles interacting only 
through gravity reads: 



where 



where 



N 



$, 



G 



N 

E 

J=0 



mj(f>{\ri 



,h). 



(3) 



(4) 



The expression for the potential $ departs from its New- 
tonian definition due to the introduction of softening and 
hence of the kernel (j> and the related length scale h. When 
this quantity has no further dependences and is kept fixed 
the resulting equation of motion will resemble the Newto- 
nian original, but for the inheritance of this additional scale; 
its expression is derived by applying the Euler-Lagrange 
equations to ^ and for particle i it will read: 



dv. 



^ -£''.('■)■ 



i=o 



where 



F,,ih) 



-VcE>, 



-GrUi 



,h)- 



(5) 



(6) 



Here cf>' = d4>/d \ri — rj\ is the force kernel. The expressions 
for both the potential and the force kernel in the cubic spline 
case can be found in BK09 (Appendix A) . The force Fij{h) 
reduces to the Newtonian gravitational force for separations 
\ri — rj\ > h and to the kernel-dependent "softened" version 
otherwise. 

When h is allowed to vary from particle to particle according 
to the density of the environment a few complications arise. 
In order to maintain its translational invariance and hence 
ensure that the resulting dynamics would conserve linear 
momentum, the Lagrangian needs to be symmetrised. One 
way to do this is by averaging over the potentials evaluated 
with the softenings of the interacting pair of particles; the 
new Lagrangian would read: 



L = 






N N 



1=0 j=0 



<f)ij{hi) -f (l)ij{hj 



-(7) 



where (j)ij{hi) = (l>{\ri — rj\ ,h). When deriving the equa- 
tion of motion from ((Tj), the additional dependence on the 
position r through h results in an extra term besides the 
classical inverse square law. The evolution of a particle i 
subject only to the action of gravity and supplied with a 
variable softening length will in fact obey: 



dvi 
Ht 




(8) 



^ _ dhi sr-^ 



ruk- 



d (j>ik{hi) 
dhi 



«.=■-£!: 



dhi v^ dWik{hi 



dpi 



TUk- 



dhi 



(9) 



(10) 



The additional contribution, which we will refer to as the 
"correction term" , is attractive in nature and will therefore 
act towards increasing the resulting gravitational force. If 
this term were not accounted for when using adaptive soft- 
ening, the particles would be evolved according to a law 
incoherent with the system's Lagrangian, i.e. energy conser- 
vation would be lost. 



3 IMPLEMENTATION IN GADGET 

We have implemented the adaptive softening formalism on 
the latest version of the GADGET code, namely GADGET- 
3. GADG ET-3 computes gravitational forces via the TreePM 
method (|Xulll995l : iBode. Ostriker fc XullJOOol : iBaglal [200^ ) 
and hydrodynamical interactions by means of SPH; it differs 
from the previous versions of the code in that it features a 
more flexible domain decomposition, something that made 
it suitable for simulations in volving extreme levels of cl us- 
tering, as the Millennium-II (jBovlan-Kolchin et al.ll2009l ). 
Implementing adaptive softening required modifications of 
the Tree- algorithm and of the timestep criterion along with 
the introduction of the machinery to compute the soften- 
ing lengths. Besides the obvious change in the expression 
for the gravitational acceleration, the opening criterion for 
the nodes has been modified in order to avoid the pres- 
ence of smoothed particle-node interaction; in other words, 
the correction term can be non-zero only within a particle- 
particle interaction (more details in Sec. 3.2.4 of BK09). The 
timestep criterion employed by GADGET for coUisionless par- 
ticles reads: 



Ai = 



2776 



1/2 



(11) 



where -q is an accuracy parameter, e is the Plummer equiva- 
lent softening {h ~ 2.8e, where h is the support of the cubic 
spline kernel) and a the acceleration of the particle. We kept 
the criterion itself unchanged and substituted e with the in- 
dividual value of the adaptive softening. 
As to the computation of the softening lengths the method is 
identical to the one GADGET uses for setting the smoothing 
lengths in SPH calculations. Qualitatively, they represent 
the radius of the sphere centred on the particle and con- 
taining a specific number of companion particles, generally 
referred to as "neighbours" ; formally, this translates into the 
following relation: 



N 
3 V^ 



^h'^'^W,j{h,) = N„gbs (12) 

i=o 

where N^gbs stands for the number of neighbours and N is 
the number of particles within distance hj from the target 
particle j. The above equation is solved iteratively with a 
Newton-Raphson method until the difference between the 
two sides falls below a certain tolerance threshold. 
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Adaptive softening lengths can be activated both when the 
code works in TreePM mode and when it uses the Tree- 
only algorithm. In the latter case the softening lengths are 
left varying without boundaries according to the local fea- 
tures of the particle distribution; in the former case we do 
instead allow for the presence of a minimum and maximum 
value for the softening lengths. As explained in Sec. 3.2.1 
of BK09, the existence of a minimum value is not crucial 
and only prevents the simulation from becoming overly ex- 
pensive in terms of computational time; conversely, the up- 
per bound is introduced to ensure that the long-range force 
(the particle-mesh contribution) is negligible on the scales 
where softening is important, so that errors arising from the 
non-modification of the long range force are under control. 
These bounds are expressed in terms of the splitting scale 
Ts , the scale (generally of order the grid spacing) where the 
splitting of the potential in a long-range and short-range 
component is performed; choosing hmax — ■^ results in the 
long-range contribution being below 1% of the total force 
at scales where softening is important. Although we always 
impose a lower limit to the softening length when using the 
code in its TreePM mode, we did not find the presence of 
an upper limit to have dramatic consequences on the re- 
sults, especially when using the adaptive formalism in its 
full, conservative version. 



4 TESTS 

In this section we present some of the tests performed 
in order to check the correctness of the implementation 
and explore the general effects of adaptive softening when 
simulating different physical scenarios. We will initially 
show the behaviour of the code in simulating simple 
systems, of well-known properties; the force profile of a 
Plummer and Hernquist models are investigated, and their 
temporal evolution in and out of equilibrium; the density 
profile of a polytrope and the behaviour of its total energy 
in time is also shown. Most of these examples are present 
already in PM07 and were specifically chosen to test our 
implementation. 

In all the numerical simulations presented in this section 
we adopt units of mass [M] = 1, length [i?] = 1 and G = 1. 
As a result, the energy per unit mass is measured in units 
of GM/R and time in [GM/ R^)-^^'^ . 



4.1 Evolution of a Plummer sphere 

The system considered here consists of a set of A'^ particles 
distributed according to a "Plummer" profile: 



p(r) = 



3GMr? 



47r(r2-f r2)5/2 



(13) 



Here the total mass M and the scale radius r^ are set equal 
to unity. The idea is to evaluate the resulting gravitational 
force profile and investigate its dependence on the choice of 
softening. Once this is accomplished we concentrate on the 
behaviour of the total energy as the system is let evolve in 
time. 

This test is identical to that presented by PM 07 and we refer 
to section 4.3 of their paper or alternatively to lAarseth et al.l 
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Figure 1. Average squared errors (ASE) of the force field gener- 
ated by distributing A'^ = 1000 particles according to a Plummer 
profile. The left panel shows the behaviour of the ASE as a func- 
tion of the choice of (fixed) softening; the right panel shows the 
results of adaptive softening varying the number of neighbours. 



llQTJ) for details on the setup of the initial conditions. 
The test was run using different number of particles and 
both Tree-only and TreePM algorithms for the evaluation 
of the gravitational force; the results behave as expected in 
the different cases and here we show only those obtained 
using the pure Tree method on A'' = 1000 particles. 
Fig. [T] shows the averaged square errors (ASE) of the simu- 
lated force field corresponding to different choices of both 
fixed and adaptive gravitational softening. This quantity 
measures the deviation of the force experienced by particles 
at different radii from the analytical value, given by 



f{r) 



GMr 



(^2 _!_ ^2)3/2 ' 

and it is defined as 



ASE: 



1 



f 2 /v A^ 



I Ji J exact \J^i) \ 



(14) 



(15) 



where fi is the force on particle i and fmax is the maximum 
value of the exact solution. For a discussion on the use 
of the ASE or related quantit ies to assess t he er ror on a 
force profile we refer to PM07, iMerrittI (119961 ) and iDehnenI 
(|200l[ ). 

The results show how less sensitive the representation 
of the force field is to the choice of Nngbs than to the 
choice of e; choosing from 30 to 500 neighbours does not 
shift the global ASE much away from the minimum value 
corresponding to the "optimal" choice of fixed softening 
(topt ~ 0.2, of the order the mean interparticle separation 
within the scale radius). As already commented by PM07, 
the slightly worse behaviour obtained when the new force 
definition is used and Nngbs < 100 can be ascribed to 
noise-induced gradients in the softening lengths altering the 
evaluation of the correction term; notwithstanding this, the 
ASEs for the "adaptive-|-correction" case are remarkably 
close to the minimum corresponding to the optimal choice 
of fixed softening throughout the full range of Nngts ■ These 
results compare very well with those of PM07 (see their 
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Figure 2. Behaviour of the total energy (kinetic plus potential) 
of the Plummer sphere as a function of time. The velocities were 
generated according to the equilibrium distribution function. For 
a description of the units, sec the introduction to Sec. |4] 



Evolution of the (expanding) Plummer sphere 



i -0.066 



1 " 


■ " 1 1 






Adoptive 




/ 
/ 
/ 

. / 

' * 

7 


- 







30 
Time 



Figure 3. Behaviour of the total energy (kinetic plus potential) 
of the Plummer sphere as a function of time. The equilibrium 
velocity distribution has been multiplied by a factor 1.2, leading 
to an initial global expansion of the system. For a description of 
the units, sec the introduction to Sec. |4] 



Fig. 2 and the second set of lines from the top). Choosing 
the optimal value for e and N^gbs = 60 we have evolved the 
system in time and checked the conservation of the total 
(kinetic plus potential) energy. Fig. [2] and [3] show the results 
for two different initial velocity setups, corresponding to a 
dynamical equilibrium and a perturbed state respectively. 
In the first case we expect no major evolution of the 
system (at least for several relaxation times) whereas in 
the second an initial overall expansion occurs before a state 
of equilibrium is re ached. We refer again to PM07 and 
lAarseth et al.l (jlOTJ) for details on the setup of the velocity 
profiles. The energy fluctuates or increases, reflecting the 
global changes in the softening lengths, when the adaptive 
formalism is used without changing the equation of motion 
accordingly; we register fluctuations of order per cent in 
the equilibrium set-up and an increase of « 6% in the 
perturbed set-up. Conservation is instead re-established 
down to time-stepping accuracy as soon as the correct 
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Figure 4. Average squared errors (ASE) of the force field gener- 
ated by distributing A'^ = 1000 particles according to a Hernquist 
profile. The left panel shows the behaviour of the ASE as a func- 
tion of the choice of (fixed) softening; the right panel shows the 
results of adaptive softening varying the number of neighbours. 



equation of motion is used. Note that the initial total 
energies are different in the adaptive and fixed softening 
cases; this is due to the definition and evaluation of the 
potential energy being different in the two scenarios (see 
Sec. Ell. 



4.2 Evolution of a Hernquist model 

We now perform a similar analysis on a different density dis- 
tribution; we generate a Mon te Carlo realisation of a Hern- 
quist model (|Hernauistll 19901 '). whose density profile reads: 



p(r) = 



GMrs 



(16) 



27rr(rs -|- r)^ 

The model is cusped near the origin and provides a more re- 
alistic representation of the matter distribution in collapsed 
objects of cosmological interest, but the steep rise in the 
density profile makes the distribution more difficult to re- 
solve than the Plummer case addressed in the previous sec- 
tion. We investigate the ASE and the energy evolution for 
a Hernquist model in equilibrium, sampled with A'^ = 1000 
particles; as in the previous section, the total mass M and 
the scale radius r^ are set equal to unity. Fig. [4] shows the 
averaged square errors (ASE) of the simulated force field 
corresponding to different choices of both fixed and adaptive 
gravitational softening. Again we see how adapting softening 
provides relatively stable results when varying the number 
of neighbours from 30 to 500; conversely, one needs to choose 
the fixed softening accurately to avoid a strong misrepresen- 
tation of the force profile. We note that the optimal value for 
e is smaller here than for the Plummer sphere {topt ~ 0.05 
vs. eopt ~ 0.2), a consequence of the different density distri- 
bution in the two models: the Hernquist being cusped and 
the Plummer being flat at the centre. A number of neigh- 
bours Nngbs ~ 60 is enough to overcome the noise in the 
evaluation of the correction term and to provide better re- 
sults than in the case where plain adaptive softening is used. 
For any number of neighbours, however, the force profile of 
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Figure 7. Cumulative density profile for the 1000-particlcs Hernquist model at four different times (i = 200,250,300,350). The gray 
lines represent the initial profile at t = 0; solid lines correspond to the fiducial choice for e(0.05) and A'^na()s(60); dotted lines correspond 
to e = 0.001 and Af„gi,g = 30, while dashed lines to e = 1 and N^gbs = 500. 



the Hernquist model is reproduced better by adapting the 
gravitational softening. The results are in qualitative agree- 
ment with those of PM07 (see their Fig. 3 and the second set 
of lines from the top). As before, we now choose the optimal 
value for e and N^gbs = 60 to evolve the system in time, 
checking the behaviour of the total energy. The results are 
shown in Fig. [5] for the equilibrium model; again, the energy 
fluctuates (« 1%) when the adaptive formalism is used with- 
out changing the equation of motion accordingly, whereas 
conservation is re-established down to time-stepping accu- 
racy as soon as the correct equation is used. 
A more demanding test for energy conservation is to follow 
the collision of two such particle distributions and check if 
the adaptive formalism maintains its conservative proper- 
ties even in this more dynamically active scenario. To this 
purpose we have simulated the head-on collision of two dif- 
ferent realisations of the Hernquist model, where the two 
systems were placed at a distance r of 150 scale radi i and 
made approach each other at a speed v = yj {2GM/r). The 
evolution of the total energy is shown in Fig. [S] the collision 
between the two spheres at t ~ 900 leaves a clear imprint 
in the total energy, when the softening is allowed to vary 
and the standard equation of motion is used, while as soon 
as the correction term is taken into account conservation is 
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Figure 5. Behaviour of the total energy (kinetic plus potential) 
of the Hernquist model as a function of time. The velocities were 
generated according to the equilibrium distribution function. For 
a description of the units, see the introduction to Sec. |4] 



recovered. 

We now want to investigate whether or not the use of adap- 
tive softening helps maintaining the original slope of the 
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Figure 6. Behaviour of tlie total energy of two colliding Hern- 
quist spheres as a function of time. The collision is head-on, with 
the two systems at an initial separation r = 150 r^ and approach- 
ing each other at a speed v = •v/(2GM/r). For a description of 
the units, see the introduction to Sec. |4] 



Figure 8. Radial density profile of a polytropic sphere in hydro- 
static equilibrium. The solid black line represents the numerical 
solution of the Lane-Emden equation (Eq.[T7]l; the points are the 
results obtained by evolving a N = 1445-particles system with 
different choice of softenings. 



density profile in time. To this purpose we let the equilib- 
rium system evolve for several dynamical times and compare 
the matter distribution among the different models. Fig. [7] 
shows the cumulative density profile at four different times, 
marked in the upper left corner. In all the panels the gray 
line represent the initial profile; the solid lines correspond 
to the optimal choice for fixed softening {eopt ~ 0.05) and 
to Nngbs ~ 60 for adaptive softening; tlie dotted lines cor- 
respond to e = 0.001 and Nngts = 30, whereas the dashed 
ones to e = 1.0 and Nngbs = 500. Again we see how sensi- 
tive the results are to the choice of e and, correspondingly, 
how little they depend on Nngbs , especially when the correct 
equation of motion is used. The density profile in the runs 
with adaptive softening is compatible with the original one 
at all times (and for all Nngbs, when the correction term is 
added), with a tendency to outperform the results obtained 
when using fixed softening. 



4.3 Equilibrium structure of a polytrope 

Here we consider a system evolving under the action of 
both gravity and hydrodynamical forces. A homogeneous 
gas sphere of initial radius ro = 1 with equation of state 
P = Kp'' , 7 = 5/3 and initial null internal energy is released 
to the action of self-gravity and pressure force until hydro- 
static equilibrium is reached. The process is speeded up by 
means of the standard SPH artificial viscosity together with 
a damping term in the force equation; their effect is to de- 
liberately remove kinetic energy helping the system settling 
faster to equilibrium. An analytical solution for the internal 
structure of the system does not exist when 7 — 5/3; the ex- 
act radial density profile was then obtained by numerically 
integrating the corresponding Lane-Emden equation: 



jK 



d' 



47rG(7 - 1) dr^ 



{rp^ 



+ rp = 0. 



(17) 



The reference for the setup and realisation of the test is Sec. 
4.4 of PM07 and there we refer the reader for additional de- 
tails on the generation of initial conditions and features of 



the run. 

Fig. |8] shows the density profile of the simulated system at 
t = 40, plotted against the equilibrium solution given by 
Eq. 1171 The units are as those introduced in the previous 
section; the time was chosen by following the oscillations in 
the density at r = 0.2: when their amplitude reduced of a 
factor of 20 with respect to the initial, we assumed equilib- 
rium was reached. 

Three cases are shown, namely fixed softening (e ~ 1/40 of 
the mean interparticle separation, as chosen by PM07) and 
adaptive softening {N^gbs ~ 60) with and without the use 
of the correction term. In all cases we have used 60 neigh- 
bours for the SPH computations. As found by PM07, the 
use of adaptive softening provides a better representation of 
the density profile, especially the inner regions; the higher 
central densities can be related to a smaller gravitational 
softening than in the fixed case (at least when the correc- 
tion term is employed), while the slightly better agreement 
at large radii comes from much larger softenings reducing 
the noise of the particle distribution. 

The equilibrium solution just found has then been per- 
turbed by inducing radial oscillations on the sphere. The 
system was let evolve under gravity and pressure force only, 
the behaviour of its total energy being now of interest. The 
results are displayed in Fig. |9l The system oscillates with a 
period r « 4, close to the expected 3.82 for the fundamental 
mode of oscillation of such a polytrope. This is seen clearly in 
the behaviour of the energy when adaptive softening is used 
and the correction to the equation of motion is neglected; the 
evaluated potential becomes cyclically deeper and shallower 
following, respectively, the overall decrease and increase of 
the softening lengths. Energy conservation is re-established 
at the level of the runs with fixed softening once the system 
is evolved according to the appropriate equations. 
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Figure 9. Behaviour of the total energy (kinetic, potential and 
internal) of the polytropic sphere as a function of time. The sys- 
tem has been perturbed from the equilibrium state at t = 40 by 
the induction of radial oscillations. For a description of the units, 
see the introduction to Sec. |4] 
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Figure 10. Number of particle with (^ (Eq. |9ll in a certain range 
as a function of f . The code was run on the same clustered distri- 
bution of particles using different values for the number of neigh- 
bours. 



5 PERFORMANCE IN A COSMOLOGICAL 
ENVIRONMENT 

The problem of the formation and evolution of the large- 
scale structure of the Universe has no analytical solution 
and this complicates the analysis when different simulations 
of the same system are compared; to assess the results and 
determine which technique is dealing best with the prob- 
lem is not obvious anymore. Increasing the resolution of a 
cosmological simulation generally goes along with extending 
the representation of the initial perturbation field to smaller 
scales, resulting in new fluctuations entering the initial con- 
ditions and therefore, strictly speaking, in a new system. In 
order to address the problem of structure formation and of 
how this is dealt with by the different softening formulations 
we have decided to perform a series of simulations sharing 
the same initial power spectrum of fluctuations, but differ- 
ing in the total number of particles (|Binnevll2004l 'l : the same 
structures form at all the resolution levels, but are more 
accurately described as more and more particles populate 
them. If the use of adaptive softening allowed to anticipate 
the behaviour shown by the standard runs at higher reso- 
lution, this would assess the superiority of the method over 
the usual choice of flxed softening. 

The simulated system here consists of a 100 h^^Mpc side 
periodic box in a ACDM cosmology defined by the following 
choice of parameters: 



fttot =1, ^n =0.3, Qi, = 0.04, f^A 

h = 0.7, 0-8 = 0.9, ns = 1, 



0.7, 



(18) 



where h and as are the values, at redshift zero, of the di- 
mensionless Hubble parameter and rms of the mass fluc- 
tuations smoothed on a scale of 8 h~ Mpc, whereas Us is 
the index of the primordial spectrum of fluctuations. Two 
sets of simulations have been performed: three runs with 
64'^, 128'' and 256'* particles using fixed softening and two 
runs with 64'^ and 128'^ particles using adaptive softening 
(with and without the correction term). The initial condi- 
tions were generated at z = 60 and differ among them only 
in the number of particles resolving the same initial fiuctu- 



ation field. The pre-initial uniform distribution, represented 
by an equally-spaced grid of particles, is applied a displace- 
ment field generated on 64^, 128^ and 256"' grids, depending 
on the resolution level. The power spectra are anyway iden- 
tical in all the three cases and they are truncated at the 
smallest frequency resolved on a 64'^ grid. The simulations 
follow only the evolution of the dark matter component of 
the density field and the resulting mass associated to a sin- 
gle particle varies from ~ 32 to 4 and 0.5 • 10^" h~^MQ 
depending on the resolution level of the run. In all the cases 
we have used the TreePM algorithm varying the grid from 
64'^ to 128'^ and 256^ , according to the resolution level. The 
features of the simulations are summarised in Table [T] 
The value of flxed softening was chosen to be 1/40 of the 
mean interparticle separation in the system, i.e. the smallest 
advisable value according to the generally accepted criterion 
that limits this scale to avoid too negative binding energies 
in close pairs. The choice of the number of neighbours for the 
adaptive runs corresponds to the minimum number ensuring 
a robust evaluation of the correction term. To this purpose 
we have monitored the changes in the two crucial quanti- 
ties entering the evaluation of the correction term, namely 
(" and Q (Eq. [9] and 1101 respectively) , when varying Nngts ■ 



The code was run on the 



snapshot of the 128 sim- 



ulation employing fixed softening with N„gbs varying from 
16 to 80 with a step AN„gbs = 8. Fig. [10] and [11] show the 
distribution of the two quantities when varying Nngts', we 
can conclude that a number Nngts = 60 is enough to ensure 
a converged evaluation of ( and Q and that the resulting 
correction term would be robustly estimated. We have also 
run the simulations using adaptive softening and no correc- 
tion to the equations of motion, using the same number of 
neighbours as the fully adaptive ones; although the results 
of the previous section suggest that the introduction of the 
correction term is crucial to the proper functioning of the 
adaptive formalism, we will anyway show the behaviour of 
these simulations for completeness. 
In the following subsections we summarise the main results. 
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Table 1. Basic features of the simulations presented in Sec.[5l 
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Figure 12. Particles with an associated density greater than a threshold. Their total number is shown in the upper-right corner. The full 
100 ft-^Mpc side box is shown and all the particles arc plotted; some are not distinguishable from their neighbours due to their proximity, 
which in some cases can reach 0.5 ft— ^kpc. The densities were computed as the position-weighted sum of the nearest neighbours, in the 
usual SPH fashion; the same number of neighbours (60) has been used for all the simulations. The results are shown at z = 0. 



5.1 Global behaviour I - Clustering 

The level of clustering in the simulations has been first as- 
sessed qualitatively by investigating the densities achieved 
throughout the box. Fig. [12] shows the particles with an as- 
sociated density greater than 105 , 5 - 105 and 10^ times the 
average density; the results are from the redshift zero snap- 
shots of the 128 and 256 simulations. The densities have 
been computed in the SPH fashion, using 60 neighbours for 
all the runs. As one can see, the runs with adaptive soften- 



ing reach higher particles densities. Noticeably, the regions 
where this enhancement of clustering is observed correspond 
to the high density regions of the 256"^ simulation; indeed, 
moving from left to right (i.e. from the fixed softening case 
to the fully adaptive one) the same areas can be seen to be 
more and more populated. The numbers in the upper-right 
corner corresponds to the total number of particles surviv- 
ing the density threshold; interestingly, the number of these 
high-density particles in the 256'^ simulation correspond to 
a mass of « 23000, 10000, 300 particles at the 128^ resolu- 
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Figure 11. Number of particle with Q, fEg. llOII in a certain range 
as a function of r2. The code was run on the same clustered distri- 
bution of particles using different values for the number of neigh- 
bours. 
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Figure 13. Two-point correlation function at 2; = for a series of 
simulations of a 100 /i~^Mpc side box in a ACDM universe; shown 
is the ratio to the results of the 256'^ simulation. The runs differ 
only in the choice of softening and in the total number of particles 
used to sample the initial perturbation field. The arrows indicate 
the equivalent-Plummer values of the gravitational softening in 
the standard runs. 



tion level. 

On a more quantitative level, we have also analysed the clus- 
tering by means of the two-point correlation function ^(r). 
This statistics represents the excess probability, compared to 
a uniform random distribution, of finding pairs of particles 
at a given spatial separation. Fig. 1131 shows the ratio of the 
correlation functions obtained in the different runs to the 
result of the 256^ simulation; the curves are shown at 2 = 
and starting from separations of 10 /i~^kpc, which corre- 
spond roughly to the resolution scale at the 256"^ level. Other 
than for some discrepancies at separation of several mega- 



parsec^J, the correlation functions of the three runs with 
fixed softening overlap almost perfectly down to 100 /i~^kpc. 
Below that scale the differences in resolution result in a pro- 
gressively larger amplitude. Since the behaviour of ^(r) at 
the small-r end is mainly determined by the distribution 
of particles within individual halos, we could conclude that 
the internal distribution of particles in a structure becomes 
denser as the resolution of the simulation is increased. The 
runs with adaptive softening produce a correlation function 
which is in full agreement with that of the "fixed" run at 
the same resolution down to roughly 100 /i~^kpc; at smaller 
separation the amplitude grows instead larger, approaching 
the results obtained at higher resolution when using fixed 
softening. This is particularly evident in the 64^-run using 
adaptive softening and the correction of the equation of mo- 
tion ( "Adapt -fcorr - 64'^", blue curve): not only is the am- 
plitude at small separation (10 - 50 ft^^kpc) higher than 
the other two 64^ runs, but the behaviour at scales between 
50 and 200 /i~^kpc is indistinguishable from the reference 
128'^-run (black, dashed curve). Both the "adaptive" runs 
at the 128'^ resolution level behave well in reproducing the 
correlation function of the 256'^ simulation, the one with cor- 
rection (purple curve) even slightly exceeding the results for 
the 256'^ case on scales below 100 /i~^kpc. 



5.2 Global behaviour II - Mass function 

The search for structures in the si mulations has been car- 
ried out with the a lgorithms fof (| Davis et al.l 1 19851 ) and 
SUBFIND ISpringel et al. 2001). Structures are first identified 
as collections oi N > Nmin particles separated by mutual 
distances smaller than some fraction b of the mean inter- 
particle separation. These so-called "FOF halos" are later 
examined by SUBFIND for the identification of self-bound 
substructures and the removal of spurious background par- 
ticles. In the simulations presented here the halos have been 
searched using a linking length b = 0.16 and a minimum 
threshold of 32 particles. The internal structure of these 
candidate objects has then been probed in order to identify 
local, gravitationally-bound overdensities; those containing 
at least 20 particles were addressed to as "subhalos" leaving 
the others as part of the smooth halo component. Finally, 
particles not gravitationally bound to any substructure of 
the parent FOF halo were dismissed. Note that SUBFIND 
was modified in the unbinding part so that the evaluation 
of the gravitational potential takes into account the individ- 
ual softening length of the particles. Fig. 1 141 shows the mass 
function of the FOF halos at 2; = 0; the number of objects 
per logarithmic mass interval is plotted for the different runs 
down to the mass limit corresponding to 32 particles. When 
the correct equation of motion is used, the agreement be- 
tween the simulations with adaptive and fixed softening at 
the same resolution level is striking throughout the full mass 

^ These discrepancies have no physical meaning and are due to a 
somewhat less accurate determination of the correlation function 
at large separations. The method wc use is Monte-Carlo based 
and tuned to provide very solid results at small separations, at 
the expense of accuracy on large scales. We checked that the 
differences we register are within the typical uncertainties due to 
the method when the correlation function is computed at these 
separations. 
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Figure 14. FOF mass function at 2 = for a series of simulations 
of a 100 /i~^Mpc side box in a ACDM universe. The upper x-axis 
displays the number of particles in the 128'^ run corresponding to 
the masses in the lower axis. The mass functions are truncated 
at a mass corresponding to 32 particle masses in every resolution 
level (i.e. at R^ lO", 1.3- lO^^^ 1.7 ■ IO^/i-^Mq for the 64^,128^ 
and 256^ run, respectively). The simulations differ only in the 
choice of softening and in the total number of particles used to 
sample the initial perturbation field. 



range; when the correction term is instead ignored, the mass 
functions drop visibly at low masses, the number of objects 
containing a number of particles less than of order Nngbs 
being significantly underestimated. Overall, the runs with 
adaptive softening and correction underestimate the total 
number of FOF halos by ~ 3%, as opposed to a ~ 25% 
for those without the correction term; if we considered only 
halos containing more than 100 particles, these percentages 
would fall to ~ 1% and ~ 14%, respectively. Similar results 
are obtained when considering the "spherical overdensity" 
masses M/Jf|; this hints to the fact that the global shape of 
the structures in the simulations is not changed significantly 
by the adoption of the adaptive softening formalism. 
The upturn in the curve for the 256'^ simulation at the low- 
mass end resembles the effect of s purious halos se e n in w arm 
dark matter simulations (see, e.g. I Wang fc White! (120071) and 
their Fig. 9); these objects form from the artificial fragmen- 
tation of filaments and have masses much smaller than the 
free-streaming mass of the model being adopted, their origin 
being entirely due to the specifics of the pre-initial condi- 
tions. Having our simulations a truncated power-spectrum, 
we might in principle be hitting the same problem and have 
a contribution to our mass functions coming from such spu- 
rious structures. This could indeed be the case for the 256'^ 
simulation at the very low-mass end, a regime we are not in- 
terested in for our comparison. As for the 128'^ simulations, 
the problem would affect objects with masses corresponding 

^ These are defined as the masses of the spherical regions centred 
on the potential minimum of the smooth halo and corresponding 
to an overdensity of A, typically ^ 200, with respect to either the 
critical density Pcrit or to the mean background density pm = 
^mPcrit- Hereafter, when referring to either the virial radius or 
the virial mass of an object, we assume the overdensity to be 
defined with respect to the mean background density. 



V 10" 




100 
(kpc/h) 



1000 



Figure 15. Cumulative mass distribution for the most massive 
FOF group at 2 = 0. 



to a handful of particles that do not survive the 32-particles 
limit and therefore do not enter the evaluation of the mass 
functions. We also specify that this problem does not in- 
validate the results we displayed in the previous section in 
terms of correlation functions; the number of particles in the 
alleged spurious structures in the 256'^ simulation is insuffi- 
cient to affect the evaluation of the correlation function at 
the scales we are interested. 



5.3 Internal halo properties 

The results in terms of correlation function and mass func- 
tion suggest that the main differences between the runs are 
likely to manifest in the internal structures of the collapsed 
objects. We will investigate here what effects the different 
definitions of softening have in the representation of the 
most massive halo formed in the simulations. Fig. llSl shows 
the cumulative mass distribution in the halo out to the virial 
radius. The values differ by fractions of per cent in the var- 
ious runs and assess around 2.4 h~^Mpc. Again it is evi- 
dent how the use of adaptive softening enhances the cluster- 
ing of particles anticipating the behaviour of higher resolu- 
tion simulations. The results of the two runs with adaptive 
softening and correction ("Adapt-fcorr - 64^", blue curve; 
"Adapt-fcorr - 128'^", purple curve) are particularly worth 
of notice: the first reproduces the higher-resolution results 
down to 20 h~^kpc, whereas the second is perfectly com- 
patible with the highest resolution run down to less than 
5 h~^kpc. Another way to interpret this result is in terms of 
mean inner density as a function of the en close d number of 
artic les. As noted bv lMoore et al.l (|l998h and lPower et al.l 
2003), obtaining robust results in regions closer to the cen- 
tre, where the density is higher, de mands increas i ngly l arge 
particle numbers. As in Fig. 14 of IPower et al.l (|2003l l we 
have compared the mean inner density measured at differ- 
ent fractions of the virial radius as a function of the enclosed 
particle number for all the runs. The results are displayed 
in Fig. 1161 The 64^ simulation using fixed softening provides 
converged results at most down to ~ 6% of the virial radius. 
At smaller radii the results start to diverge from those ob- 
tained at higher resolution and cannot be trusted anymore. 
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Figure 16. Mean inner density as a function of the enclosed 
number of particle at different fractions of the virial radius. The 
halo under consideration is the same as in Fig. 1151 The solid lines 
separate the regions, at their right, where the average coUisional 
relaxation time (trez) exceeds some fraction of the age of Universe 
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Figure 17. Average subhalo mass function at z = for the 
five most massive halos in the simulations. The subhalo masses 
have been normalised to the virial mass of the host and vary 
from 8 • 10^^ to 4 ■ W^^h-^MQ, corresponding to » 20 and 1000 
particles, respectively. 



The two runs with adaptive softening behave instead very 
weU down to 3% of the virial radius. Similarly, the 128^ sim- 
ulation is reliable down to ~ 2% of the virial radius when 
using fixed softening and down to 1% of the virial radius 
when adaptive softening is employed. IPower et al.l (|2003l ) 
relate the converged side of the enclosed-p vs. enclosed-N 
plot to the regions where the average coUisional relaxation 
time trei exceeds some fraction of the age of Universe to 
(between 0.6 and 1). In our case, depending on whether we 
consider the results at 6% of the virial radius converged or 
not, we could extend the regions down to the radii where 
trei ~ 0.4to. In any case, it is not clear whether the 256'' run 
can be considered converged according to this criterion. It 
is anyway worth of notice that the 128'' adaptive run with 
correction reproduces perfectly the enclosed density of the 
256^ at 1% of the virial radius. 

Similar results hold when investigating the properties of sub- 
structures. Fig. 1 171 shows the number of subhalos with mass 
greater than a certain value; the curves represent the aver- 
age over the five biggest halos in the simulations and the 
masses have all been normalised to the virial mass of the 
host. Due to the poverty of substructures we are not show- 
ing the results for the lowest resolution simulations and we 
are in general limited to qualitative considerations. As can 
immediately be seen though, the runs with adaptive soften- 
ing lie closer to the higher resolution simulation than does 
the run with fixed softening. 



5.4 Comments 

All the results shown so far hold at redshift zero. As shown 
in Fig. 1191 even at this redshift only « 1% of the particles 
have softening smaller than 1/40 of the mean interparticle 
separation and at higher redshifts this behaviour becomes 
even more extreme. This has the effect of reducing the am- 
plitude of the correlation function more rapidly than in the 
standard runs when going back in time. A more quantita- 
tive representation of this effect will be given in the next 



section (Fig. I20|l . The mass functions are not substantially 
affected though, at least not for objects more massive than 
the 32-particles limit (Fig.EU- 

As the particle timesteps depends on the value of the grav- 
itational softening (see Eq. Ill|l , it is natural to expect that 
they would now span a wider range of values. Fig. 1181 shows 
that this is the case; plotted are the distributions of particles 
in time binqj at redshift zero for the three 128"' simulations: 
the "adaptive" runs tend to have particles in lower bins than 
in the "fixed" simulation and, at the same time, more parti- 
cles in the highest bin. This is a consequence of what was just 
mentioned, namely that only a very small fraction of the par- 
ticles obtains an adaptive softening smaller than the fiducial 
choice for the standard simulations at the same resolution 
level. As finer time bins get populated, the overall number of 
timesteps increases when adaptive softening is used (along 
with the number of operations performed within them); at 
the same time, since the highest bins become more crowded, 
the average number of active particles - i.e. those that are 
advanced in a timestep - decreases instead, almost compen- 
sating for this overhead. Increasing Nngba has the effect that 
one intuitively imagines, namely that of increasing the soft- 
ening associated to the particles; the differences again stand 
out more evident when looking at the correlation functions, 
whose amplitude at small separations reduces, whereas the 
mass functions are confirmed to be considerably less sensi- 
tive to variations in Nngbs- 
All the simulations discussed in this section the gravita- 
tional softening was prevented to fall below 0.1% of the 
TreePM splitting scale Ts (corresponding to 558 pc for the 
64'' case and to 279 pc for the 128'' one; in both cases, this 



^ In GADGET the simulated timespan is mapped onto the inte- 
ger interval [0,2^*"'"'='''"=]. This interval is split recursively into 
timebins, the smallest of which has length 2. Each time, the code 
computes individual timesteps for the particles and distributes 
them in the correspondent bin; the smallest populated bin sets 
the next timestep for the simulation. 
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Figure 18. Cumulative distribution of particles in time bins at 
z = 0. The simulations wore run using 29 time bins. Smaller 
numbers corresponds to finer intervals. 
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Figure 19. Cumulative distribution of softenings at 2 = 0. The 
softenings are expressed in terms of equivalent-Plummer values 
(e). The red (blue) dashed curve represents the value of the grav- 
itational softening in the 128 (64 ) "fixed" simulations. 



translates to « 1.5% of the softening scale in the "fixed" 
runs); this implies that the maximum density achievable in 
the simulation corresponds to roughly 10^ p. As pointed out 
in Sec.[3l the presence of a lower bound is not crucial and it 
is introduced mainly to prevent few particles in over-dense 
regions from slowing down the simulation. No upper limit to 
the gravitational softening was set in these simulations; we 
do not find substantial differences between the runs with and 
without such a constraint, especially if the correct equation 
of motion is used. We have also checked both sets against 
Tree-only runs and found perfectly compatible results. 
Although we both observe an increased level of clustering in 
massive objects, our global results in terms of particle corre- 
lation function and halo mass function somewhat differ from 
those of BK09; when comparing the results from the adap- 
tive runs to the standard case with fixed softening, BK09 
register a deficit in the number of small-mass objects and a 
slight lowering in the amplitude of the two-point correlation 
function at small separation, whereas we notice no change 



in the mass function and instead an overall enhanced level 
of particle clustering. The simulated background cosmology 
differ, but this should have no effect in the analysis we are 
interested here, which is focused on relative differences be- 
tween the runs more than on absolute results of cosmologi- 
cal interest. The larger linking length we have employed is 
not at the origin of the discrepancy between the FOF mass 
functions, nor we think is the different evaluation of the cor- 
relation function responsible for the antithetical outcomes. 
The results differ also on the more technical timing level; 
in our simulations we do not register substantial saving of 
computing time when adaptive softening is used, but we do 
notice increasingly better performances when the number 
of neighbours is increased. BK09 halve the wallclock time 
of the reference run when using adaptive softening with 32 
neighbours and leave it unchanged when using 48. If we set 
an upper limit to the softening of order the splitting scale 
Ts (as in BK09), we progressively cut the timing when going 
from 32 up to 60 neighbours and leave it essentially unal- 
tered by increasing the number even further; the minimum 
wallclock time is reached for N^gts — 60 and it equals that 
of the simulation with fixed softening. If no upper limit is 
set for the softening, the wallclock times tend to increase; a 
minimum is reached for a number of neighbours ~ 40 and it 
is around ~ 30% higher than in the standard run. 
We ascribe the responsibility for the different behaviour of 
the implementation in the various aspects to the respective 
mother codes. 



6 SIMULATING A MINI VERSION OF THE 
MILLENNIUM-II 

As a more demanding application, we have investigated 
the effect of adaptive softening on a cosmological simu- 
lation whose small-scale clustering is known in detail by 
an extremely high resolution run. We have used the mod- 
ified GADGET-3 on the initial conditions of the "mini"- 
Millennium-II simulation (hereafter mmll), a low-resolution 
version of the better know n Millennium-II (hereafter mil, 
iBovlan-Kolchin et all [20091 ). The mmll follows the evolu- 
tion of 432'^ particles (as opposed to 2160'^ particles for the 
mil simulation) within a box of side 100 h^^Mpc. The un- 
derlying cosmology is a ACDM with parameters 



h-- 



t =1, t7„ =0.25, Qb = 0.045, Qa = 0.75, 
: 0.73, as = 0.9, n^ = 1 . 



(19) 



This results in the particles having masses m = 8.61 ■ 
10^" h~^MQ , a factor 125 larger than in the mil simula- 
tion, corresponding to the s ame mass resolution of the orig- 
inal Millennium simulation (jSpringel et al.ll2005l ). 
The initial conditions for the mmll are identical to those for 
the main (mil) run. This means that although they share 
the power spectrum of the fiuctuation field, this is sampled 
with a factor of 125 less particles in the case of the mmll; 
it is reasonable to expect that the smallest perturbations 
are then not represented. This reverses the situation with 
respect to the series of simulations analysed so far in this 
section. Fig. 1201 shows the ratio of the correlation functions 
from the mmll and "adaptive" -mmll to the original mil at 
four different redshifts. At high redshifts, down to approxi- 
mately z — 2, the amplitude of ^(r) at scales smaller than 
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Figure 20. Comparison of the two-point correlation functions 
from the mil (black curve), mmll (green curve) and "adaptive"- 
mmll (purple curve) simulations at different redshifts; shown is 
the ratio to the results from the mil. See the text for a description 
of the three cases. 



Figure 21. POP mass functions from the mil (black curve, dia- 
monds), mmll (green curve) and "adaptive" -mmll (purple curve) 
simulations at different redshifts. See the text for a description of 
the three cases. 



~ 20 h~^kpc is lower in the "adaptive" -mmll than in the 
mmll; the result is not surprising though, as at those red- 
shifts we expect almost all the particles to have softenings 
considerably larger than e = 5 h^^kpc, the value adopted 
in the original mmll. At low redshifts the situation reverses, 
leading to a correlation functions which approaches more 
and more tightly the results of the mil simulation; at z = 
and on scales below ~ 10 h^^kpc, the "adaptive" correlation 
function is a factor ~ 1.3 above that of the mmll. Fig. [51] 
shows the mass functions for the three simulations at the 
same redshifts as in Fig. 1201 The halos are identified by the 
FOF algorithm, using 6 — 0.2 and a minimum of 20 particles. 
Again, the mass function are in agreement at all redshifts. 
As for the tests presented before, this demonstrate that us- 
ing adaptive gravitational softening while incorporating the 
correction term into the equation of motion allows to resolve 
the smallest scales better than simulations with fixed soft- 
ening at the same resolution. The results converge to what 
is expected from much higher resolution simulations and no 
substantial degrading of the results is observed at high red- 
shifts, where the formal resolution of the simulations using 
adaptive gravitational softening is naturally smaller than in 
the simulations with fixed softening. 



7 CONCLUSIONS 

We have implemented adaptive gravitational softening in the 
cosmological TreeP M code GADGET-3. The f ormalism was 
introduced first by IPrice fc MonaghanI (|2007l ) and features 
the same technique used by SPH to determine individual 
softening lengths for each of the simulation particles. The 
spatial variation of this scale requires a modification of the 
equation of motion governing the evolution of the particles' 
trajectories in order to be consistent with the new depen- 
dencies introduced in the system's Lagrangian. 
We have applied this technique to several test cases in order 
to check the behaviour of the total energy when the new 
equation of motion is used; we then moved to the cosmo- 
logical scenario and, specifically, to the simulation of the 



large-scale structure of the Universe, where the evaluation 
of the effects of adaptive softening is complicated by the lack 
of an analytical solution to the problem. 
Our main conclusions are: 

• The inclusion of the correction term in the equation of 
motion is essential to ensure the conservation of total energy. 

• In cosmological simulations of the large-scale structure 
a number of neighbours ~ 60 is needed to obtain a converged 
estimation of the correction term. 

• With such a choice we show that - in contrary to previ- 
ous claims in literature - the adoption of adaptive softening 
does not lead to an under-representation of halos at low 
masses, if the correct equation of motion is used. 

• Using the adaptive scheme effectively increases the dy- 
namical range of cosmological simulations while the com- 
putational costs only mildly increase. Especially the ampli- 
tude of the two-point correlation function at small scales 
and the subhalo mass function improve compared to simu- 
lations with the same number of particles and fixed gravita- 
tional softening, anticipating the results obtained in higher- 
resolution simulations. 

• The convergence of the inner density profile for the most 
massive object found in the simulations improves signifi- 
cantly when the adaptive softening and the correct equation 
of motion is used. 

• When re-simulating a low-resolution version o f the 
Millennium-II simulation (jBovlan-Kolchin et al.l 12003 ) and 
comparing the results obtained with fixed and adaptive 
softening, we notice again perfect agreement in the mass 
functions at all times and an evolution of the "adaptive" 
correlation function towards the higher amplitude of the 
Millennium-II's at late times. 

We have so far applied the method to scenarios with equal- 
mass particle species, either gas-only or dark-matter-only; 
simulations including both species with different masses 
would also considerably benefit from the adoption of the 
scheme: the adaptive behaviour of the resolution scale would 
allow to follow the collapse of dark matter and particularly 
gas down to scales currently unachievable in standard simu- 
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lations at comparable resolution. However, having mixtures 
of particles with different masses is a non-trivial extension 
of such scheme and will be investigated in future work. Such 
future studies will be of particular importance considering 
that cold gas and star particles in hydrodynamical simula- 
tions are likely to clump at lengthscales below those typical 
chosen for the gravitational softening. 
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